function [yOut_cu] = modelFit(x_cu,model)
mats   = model.maturities;
if model.pruningOn == 1
    if model.appOrder == 1
        xf_cu      = x_cu(1:model.nx,1);
        y_cu       = model.g0  + model.gx*xf_cu;
        yields_cu  = model.by0(mats) + model.byx(mats,:)*xf_cu;
        if model.inclSurveys == 1
            ER_cu  = model.r0 + model.rx*xf_cu;
        end
        x_cu = xf_cu;
    elseif model.appOrder== 2
        xf_cu                 = x_cu(1:model.nx,1);
        xs_cu                 = zeros(model.nx,1);
        xs_cu(1:model.nx1,1)  = x_cu(model.nx+1:model.nx+model.nx1,1);
        xfxf_cu               = kron3(xf_cu,xf_cu);
        y_cu                  = model.g0        + model.gx*(xf_cu+xs_cu) + model.Gxx*xfxf_cu + 0.5*model.gss;
        yields_cu             = model.by0(mats) + model.byx(mats,:)*(xf_cu+xs_cu) + model.Byxx(mats,:)*xfxf_cu + 0.5*model.byss(mats,:);
        if model.inclSurveys == 1
            ER_cu  = model.r0 + model.rx*(xf_cu+xs_cu) + model.Rxx*xfxf_cu + 0.5*model.rss;
        end
        x_cu = xf_cu+xs_cu;
    elseif model.appOrder == 3
        xf_cu                 = x_cu(1:model.nx,1);
        xs_cu                 = zeros(model.nx,1);
        xs_cu(1:model.nx1,1)  = x_cu(model.nx+1:model.nx+model.nx1,1);
        xrd_cu                = zeros(model.nx,1);
        xrd_cu(1:model.nx1,1) = x_cu(model.nx+model.nx1+1:model.nx+2*model.nx1,1);
        xfxf_cu               = kron3(xf_cu,xf_cu);
        xfxfxf_cu             = kron(xf_cu,xfxf_cu);
        xfxs_cu               = kron3(xf_cu,xs_cu);
        y_cu                  = model.g0 + model.gx*(xf_cu+xs_cu+xrd_cu) ...
                                + model.Gxx*(xfxf_cu+2*xfxs_cu) + model.Gxxx*xfxfxf_cu ...
                                + 0.5*model.gss + 3/6*model.gssx*xf_cu;
        yields_cu             = model.by0(mats) + model.byx(mats,:)*(xf_cu+xs_cu+xrd_cu) ...
                                + model.Byxx(mats,:)*(xfxf_cu+2*xfxs_cu) ...
                                + model.Byxxx(mats,:)*xfxfxf_cu ...
                                + 0.5*model.byss(mats,:) + 3/6*model.byssx(mats,:)*xf_cu;
        if model.inclSurveys == 1
            ER_cu  = model.r0 + model.rx*(xf_cu+xs_cu+xrd_cu) ...
                + model.Rxx*(xfxf_cu+2*xfxs_cu) + model.Rxxx*xfxfxf_cu ...
                + 0.5*model.rss + 3/6*model.rssx*xf_cu;
        end
        x_cu = xf_cu+xs_cu+xrd_cu;
    end
else
    if model.appOrder == 1
        y_cu       = model.g0(mvars) + model.gx(mvars,:)*x_cu;
        yields_cu  = model.by0(mats) + model.byx(mats,:)*x_cu;
        if model.inclSurveys == 1
            ER_cu  = model.r0 + model.rx*x_cu;
        end
    elseif model.appOrder == 2
        xx_cu      = kron3(x_cu,x_cu);
        y_cu       = model.g0 + model.gx*x_cu + model.Gxx*xx_cu + 0.5*model.gss;
        yields_cu  = model.by0(mats) + model.byx(mats,:)*x_cu + model.byxx(mats,:)*xx_cu + 0.5*model.byss(mats,:);
        if model.inclSurveys == 1
            ER_cu  = model.r0 + model.rx*x_cu + model.Rxx*xx_cu + 0.5*model.rss;
        end
    elseif model.appOrder == 3
        xx_cu      = kron3(x_cu,x_cu);
        xxx_cu     = kron3(xx_cu,x_cu);
        y_cu       = model.g0 + model.gx*x_cu + model.Gxx*xx_cu + model.Gxxx*xxx_cu + ...
                     0.5*model.gss + 3/6*model.gssx*x_cu;
        yields_cu  = model.by0(mats) + model.byx(mats,:)*x_cu + model.byxx(mats,:)*xx_cu ...
                    + model.byxxx(mats,:)*xxx_cu ...
                    + 0.5*model.byss(mats,:) + 3/6*model.byssx(mats,:)*x_cu;
        if model.inclSurveys == 1
            ER_cu  = model.r0 + model.rx*x_cu + model.Rxx*xx_cu + model.Rxxx*xx_cu ...
                    + 0.5*model.rss + 3/6*model.rssx*x_cu;
        end
    end
end

% Matching output to the data 
pos_l     = find(strcmp({'$l_t$'},model.labely));
pos_w     = find(strcmp({'$w_t$'},model.labely));
%pos_c_ba1 = find(strcmp({'$c_{t-1}$'},model.labelx));
%pos_myz   = find(strcmp({'$\mu _{z,t}$'},model.labelx));
%pos_c     = find(strcmp({'$c_t$'},model.labely));
%dc        = 4*(y_cu(pos_c,1)-(x_cu(pos_c_ba1,1)+model.g0(pos_c_ba1,1))+x_cu(pos_myz,1)+log(model.params.MUZss));
pos_pai   = find(strcmp({'$\pi_t$'},model.labely));
pos_dc    = find(strcmp({'$dc_t$'},model.labely));
hoursDev  = y_cu(pos_l,1)-model.g0(pos_l,1);
wageDev   = y_cu(pos_w,1)-model.g0(pos_w,1);
dc        = 4*y_cu(pos_dc,1);
infl      = 4*y_cu(pos_pai,1);
yields    = 4*yields_cu;

if model.inclSurveys == 1
    % Conditional expectations of short rates included
    yOut_cu = [hoursDev wageDev dc infl yields' 4*ER_cu']';
else
    yOut_cu = [hoursDev wageDev dc infl yields' ]';
end
%yExtra_cu = [y_cu;yields_cu];

end
